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Abstract.  Stochastic  models  are  valuable  and  sometimes  essential  tools  for  investigating  the 
behavior  of  complex  phenomena.  In  seismology,  stochastic  models  can  be  used  to  describe 
velocity  heterogeneities  that  are  too  small  or  too  numerous  to  be  described  deterministically. 
Where  analytic  approaches  are  often  infeasible,  synthetic  realizations  of  such  models  can  be  used 
in  conjunction  with  finite  difference  algorithms  to  systematically  investigate  the  response  of  the 
seismic  wavefield  to  complex  heterogeneity.  This  paper  represents  a  continuing  effort  at 
formulating  a  complete  and  robust  stochastic  model  of  lithologic  heterogeneity  within  the  crust, 
and  the  means  of  generating  synthetic  realizations;  "complete"  implies  that  the  model  is  flexible 
enough  to  describe  all  types  of  random  heterogeneity  within  the  crust,  while  "robust"  implies 
sufficiently  constrained  parameterization  that  an  inversion  problem  may  be  well-posed.  We  use 
as  a  basis  for  investigation  geologic  maps  of  crustal  exposures  and  petrophysically  inferred 
velocities.  Earlier  efforts  at  stochastic  modeling  have  focused  on  characterization  of  the 
univariate  probability  density  function,  which  is  typically  modal  (i.e.,  binary,  ternary,  etc.),  and 
the  covariance  function,  which  is  typically  fit  with  a  von  Karman  function.  Here  we  provide  a 
means  of  characterizing  the  property  of  "sinuous  connectivity"  and  for  generating  realizations 
that  possess  this  property.  Sinuous  connectivity  is  the  tendency  for  individual  lithologic  units  to 
be  continuous  over  long  and  highly  contorted  paths;  there  is  no  means  in  the  earlier  modeling  of 
either  characterizing  or  synthesizing  this  property.  We  first  demonstrate  that  binary  sinuously 
connective  realizations  can  be  generated  by  mapping  alternating  contour  sets  from  a  Gaussian- 
distributed  surface  (a  "normal  equivalent  field")  into  the  two  values  comprising  the  binary 
probability  density.  There  is  tremendous  non-uniqueness  in  this  operation,  with  wide  classes  of 
mapping  functions  and  normal  equivalent  statistics  resulting  in  model  fields  that  are  statistically 
identical.  We  infer  from  these  observations  that  the  property  of  sinuous  connectivity  can  be 
represented  by  a  simple  binary  yes-or-no  parameter. 
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I.  Introduction 

Investigations  of  the  wavefield  response  to 
heterogeneous  media  (finite  difference  or  otherwise) 
can  be  of  a  deterministic  (i.e.,  the  exact  particulars)  or 
stochastic  (i.e.,  ensemble  properties)  nature.  The 
choice  of  one  or  the  other  is  a  matter  of  scale, 
numbers,  and  resolution.  The  deterministic  approach 
is  typically  used  when  the  structures  within  the 
velocity  field  are  few  in  number  and  large  compared 
to  the  seismic  wavelength.  In  such  cases  a  well- 
posed  inverse  problem  may  be  formulated  to  estimate 
the  velocity  field  by,  for  example,  matching  a 
synthetic  wavefield  to  the  observed  wavefield  [e.g., 
Jervis  et  al.,  1995].  The  stochastic  approach  is 
necessary  when  structures  within  the  velocity  field 
are  numerous  and/or  small  compared  to  the 
wavelength.  Large  numbers  of  scattering  features 
tend  to  make  the  wavefield  complex  (i.e., 
deterministically  unpredictable)  through  both  single 
and  multiple  forward  and  backward  scattering.  Small 
scales  make  deterministic  resolution  of  structures 
difficult  or  impossible.  Recent  studies  have  used 
finite  difference  algorithms  to  investigate  the  seismic 
wavefield  response  to  synthetic  realizations  of 
stochastic  velocity  models  [e.g.,  Frankel  and 
Clayton,  1986;  Fisk  etai,  1992;  Holliger  etal.,  1993; 
Levander  et  al.,  1994a,b],  So  far  these  efforts  have 
only  been  directed  toward  the  forward  problem, 
where  a  velocity  model  is  assumed  and  the  wavefield 
response  computed.  Qualitative  comparisons  of 
seismic  data  with  finite-difference  seismograms 
generated  in  a  variety  of  realistic  circumstances 
suggest  that  the  seismogram  is  sensitive  to  variations 
in  the  stochastic  character  of  the  medium.  Efforts  are 
now  being  directed  towards  establishing 
quantitatively  meaningful  model/data  comparisons  to 
pose  the  inverse  problem  of  estimating  stochastic 
properties  of  a  velocity  field  from  observations  of  the 
seismic  wavefield. 

This  report  is  one  in  a  series  of  efforts  to  establish 
a  "robust"  and  "complete"  stochastic  model  for 
crustal  heterogeneity  through  analysis  of  mapped 
crustal  exposures.  Though  modified  by  the  act  of 
exhumation,  crustal  exposures  represent  our  only 
direct  multidimensional  sampling  of  crustal  rocks 
from  depth.  They  provide  invaluable  information  on 
lithologies  and  their  spatial  relationships.  The 
principal  goal  of  working  with  this  data  is  not  to 
establish  precisely  the  stochastic  nature  of  the  crust, 
but  rather  to  ascertain  the  types  of  stochastic  models 
that  are  appropriate.  By  "robust"  we  intend  a  model 
with  few-enough  parameters  that  an  inversion 
problem  may  be  well-posed.  By  "complete"  we 
intend  a  model  that  can  reproduce,  through  synthetic 
realization,  all  the  important  physical  properties  of 
the  field.  There  is  a  delicate  balance  between  these 
two  concerns.  If  we  have  too  few  parameters  in  the 
desire  for  robustness,  then  the  model  may  not  be 
flexible  enough  to  characterize  the  variety  of 
stochastic  morphology  observed.  If  we  favor 
completeness  we  may  require  more  parameters  than 
we  can  ever  hope  to  solve  for  in  an  inverse  problem. 


Previous  work  on  mapped  crustal  exposures  has 
focused  on  combined  modeling  of  the  covariance 
function  and  probability  density  function  (PDF)  of 
the  velocity  field  (as  inferred  from  petrophysical 
conversion  of  the  lithographic  field).  Tlie  covariance 
function  has  often  been  successfully  modeled  using 
the  von  Karman  [1948]  model  [e.g.,  Holliger  et  al., 
1993;  Levander  et  al.,  1994].  The  PDF  has  often 
been  modeled  by  a  simple  modal  distribution  (e.g., 
binary,  ternary,  etc.)  [Holliger  et  al.,  1993;  Levander 
et  al.,  1994;  Goff  et  al.,  1994],  which  reflects  the 
observation  that  mapped  crustal  exposures,  even 
morphologically  complex  ones,  typically  consist  of  a 
small  number  of  distinct  lithologic  units  of  relatively 
constant  velocity.  Synthetic  realizations  generated 
from  these  models  successfully  reproduce  many  of 
the  important  physical  properties  of  the  observed 
field,  including  the  modal  PDF,  characteristic  scales, 
structural  anisotropy,  and  fractal  dimension. 
However,  one  property  clearly  seen  on  many  of  the 
crustal  exposures  has  thus  far  remained  beyond  the 
grasp  of  stochastic  characterization;  "sinuous 
connectivity",  the  subject  of  this  report. 

Sinuous  connectivity  is  a  difficult  property  to 
define  in  words,  but  simple  to  demonstrate.  Figure  la 
is  a  digitized  section  of  the  Lewisian  gneissic  terrain 
(Scotland),  an  exposed  section  of  the  middle  crust 
consisting  primarily  of  amphibolite  dikes  (black)  and 
gneiss  (white)  (gray  =  no  data).  The  gneiss  has  a 
petrophysically  inferred  p-wave  velocity  of  6.2 
km/sec,  and  the  amphibolite  6.75  km/sec.  Levander 
et  al.,  [1994a]  formulated  a  stochastic  model  based 
on  analysis  of  the  PDF  and  covariance  structure  for 
the  grid  shown  in  Figure  1,  and  a  synthetic  realization 
from  a  refined  version  of  this  model  is  presented  in 
Figure  lb.  The  grids  shown  in  Figures  la  and  lb 
have  essentially  identical  PDF  and  covariance,  and 
those  physical  attributes  characterized  by  those 
functions  are  well-matched;  i.e.,  the  percentage  of 
gneiss  and  amphibolite,  the  characteristic  scales,  and 
the  fractal  dimension  (the  overall  degree  of 
roughness).  Nevertheless  the  comparison  is  not 
satisfactory.  The  amphibolite  in  the  geologic  map  is 
highly  connected  and  sinuous,  whereas  in  the 
synthetic  field  it  is  disconnected  and  blob-like. 

In  this  paper  we  present  a  method  for 
incorporating  sinuous  connectivity  into  stochastic 
models  for  lithologic  (i.e.,  modal)  heterogeneity,  and 
the  means  by  which  to  generate  a  sinuously 
connective  synthetic  realization  that  honors  the 
covariance  and  PDF  structure  specified  by  the 
stochastic  model.  The  success  of  this  approach  is 
demonstrated  in  Figure  lc,  which  displays  a 
sinuously  connective  synthetic  field  with  identical 
PDF  and  covariance  structure  to  the  synthetic  shown 
in  Figure  lb.  The  comparison  between  this  synthetic 
to  the  Lewisian  data  in  Figure  la  is  clearly  superior; 
we  believe  that  this  stochastic  model,  incorporating 
PDF  (modal  or  otherwise),  covariance,  and  sinuous 
connectivity  characterization,  represents  as  complete 
a  model  as  might  be  necessary  for  characterizing 
seismic  velocity  heterogeneity  within  the  earth.  The 
robustness  of  this  model,  as  applied  to  seismic 
inversion,  will  be  the  topic  of  future  investigation. 
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Figure  1.  (a)  Digitized  map  of  the  Lewisian  gneiss  complex,  Scotland  [Levander  et  al.,  1994],  Also  mapped  but  not 
shown  are  intermediate  schist,  which  comprise  only  1%  of  the  total  and  so  were  ignored  as  insubstantial.  Map  is 
78%  gneiss  (6200  m/sec  inferred  p-wave  velocity)  and  22%  amphibolite  (6750  m/sec).  The  grid  spacing  is  0.0268 
km,  with  100  rows  and  109  columns,  (b)  Synthetic  realization  based  on  PDF  and  covariance  modeling  of  the 
Lewisian  complex,  (c)  Synthetic  realization  based  on  PDF  and  covariance  modeling  of  the  Lewisian  complex  and 
including  the  property  of  sinuous  connectivity. 


628 


II.  PDF  and  Covariance  Modeling  and  Synthetics 

Details  regarding  PDF  and  covariance  modeling 
based  on  crustal  exposure  maps  are  given  in  Holliger 
et  al.  [1993]  and  Levander  et  al.  [1994a].  Heretofore 
covariance  modeling  has  been  accomplished  by 
forward-model  fitting  of  the  von  Kdrmdn  function  to 
the  data  covariance.  Here  we  apply  the  inversion 
algorithm  of  Goff  and  Jordan  [1988]  converted  to  use 
with  gridded  data.  The  method  for  generating 
synthetic  realizations  from  such  models  is  given  in 
Goff  et  al.  [1994],  In  this  section  we  provide  a  brief 
overview  of  the  salient  points  from  these  references. 

PDF  Modeling 

In  the  case  of  a  modal  field,  modeling  the  PDF  is 
robust  and  straightforward:  the  PDF  simply  describes 
the  percentage  of  each  unit  present  in  the  field.  For 
example,  in  the  Lewisian  section  described  above, 
78%  of  the  map  is  gneiss,  and  the  remaining  22%  is 
amphibolite  (Figure  2a). 

Covariance  Modeling 

The  covariance  function,  or  its  Fourier  equivalent 
the  power  spectrum,  represents  our  primary  tool  for 
characterizing  spatial  roughness  properties.  It  is 
defined  by: 

Caa(x)  =  £[/i(x,)/i(x1+x)],  (1) 

where  /i(xj)  is  a  zero-mean,  homogeneous  random 
field  specified  at  vector  location  xj,  and  £[•] 
represents  the  expectation  function  [e.g..  Feller, 
1971],  The  variable  x  is  defined  as  the  lag  vector. 
Where  h  is  sampled  on  a  m  by  n  grid  specified  by 
discrete  (column,row)  locations  (ij),  the  discrete 
covariance  function  Chh(k,l)  can  be  estimated  by 

n-k  m-l 
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Covariance  modeling  is  typically  accomplished  by 
fitting  a  parameterized  functional  model  to  the 
sample  covariance  estimated  from  the  digital 
exposure  map  using  forward  modeling  [e.g.,  Holliger 
et  al.,  1993;  Levander  et  al.,  1994a;  Holliger  and 
Levander,  1994]  or  a  least-squares  inversion  [Goff 
and  Jordan,  1988].  The  von  Kdrmdn  [1948] 
covariance  model  [e.g.,  Wu  and  Aki,  1985;  Frankel 
and  Clayton,  1986;  Fisk  et  al.,  1992;  Holliger  et  al., 
1993;  Levander  et  al.,  1994a,  Holliger  and  Levander, 
1994]  represents  a  class  of  monotonically  decaying 
(i.e.,  aperiodic)  functions,  and  includes  as  a  subset  the 
exponential  form.  The  singular  advantage  of  the  von 
Kirmctn  model  is  that  it  explicitly  includes  the  fractal 
dimension  as  a  variable.  Goff  and  Jordan  [1988] 
modified  the  von  Kirmdn  covariance  function  to 
account  for  structural  anisotropy  in  2-dimensions 
(easily  expanded  to  3-dimensions).  The  following 
parameters  specify  the  2-D  anisotropic  von  K£rm£n 
model: 

1 .  The  rms  velocity  H  is  the  average  variation 
about  the  mean  velocity. 

2.  The  lineament  orientation  8S  is  strike  of 
direction  of  maximum  correlation;  structures  will 
tend  to  be  oriented  along  this  direction. 

3.  The  scale  parameter  controls  the  rate  of  decay 
of  the  covariance.  In  one  dimension  the  scale 
parameter  is  specified  by  kQ.  In  two  dimensions  there 
are  two  principal  scale  parameters:  kn  in  the  normal- 
to-strike  direction,  and  ks  in  the  along-strike 
direction;  kn  >  ks. 

4.  The  aspect  ratio  a  is  the  characteristic  planar 
shape  of  structures,  defined  by  the  ratio  knlks. 

5.  The  Hausdorff  (fractal)  dimension  D  is  a 
measure  of  roughness. 

An  inversion  procedure  for  parameter  estimation 
also  provides  estimates  of  uncertainties.  Estimated 
Lewisian  model  parameters  are,  with  1  -a 
uncertainties:  H  =  234  ±  4  m/sec,  8S  =  85.0°  ±  0.4°, 
kn  =  24.5  ±  1.0  km’1,  its  =  6.0  ±0.3  km*1,  andD  = 
2.54  ±  0.04.  Figure  2b  displays  the  best  fit  von 
Kirmdn  model  to  the  1-D  sample  covariance  in  the 
row  direction  (Qj|j(i,0)). 
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Figure  2.  (a)  PDF  for  the  Lewisian  stochastic  model,  (b)  Comparison  of  Lewisian  data  (solid)  and  best-fit  model 
(dashed)  1-D  covariances  for  the  row  direction  (C/,j,(i',0)). 
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Synthetic  Fields 

Realizations  of  the  stochastic  model  described 
above  can  be  synthesized  by  first  generating  a 
Gaussian-distributed  field,  which  we  term  a  "normal 
equivalent"  field,  and  then  mapping  it  into  a  new 
field,  the  "model"  field,  that  conforms  to  the  PDF 
specified  in  the  model.  This  algorithm  includes  the 
following  steps  [Goff  etal.,  1994]: 

1 .  A  mapping  function  is  created  which  specifies 
the  conversion  of  a  Gaussian  PDF  with  zero  mean 
and  unit  variance  to  the  PDF  p(h)  specified  by  the 
stochastic  model.  Specifically,  for  each  possible 
value  g  sampled  from  a  Gaussian  PDF  pG(g)<  we 
find  a  new  value  h  such  that: 

*  g 

J p(ti  )dh!  =  jpG(g'W-  (3) 


For  a  simple  binary  field,  where  unit  1  has  total 
probability  p\,  and  unit  2  has  total  probability  p2  =  1 
-  pi,  equation  (3)  is  equivalent  to  identifying  a  cutoff 
Gaussian  value  gc,  where  every  value  less  than  gc  is 
mapped  to  unit  1,  and  every  value  greater  than  gc  is 
mapped  to  unit  2,  and 


P\-  I  Pc(g'W  ■  (4) 


For  the  Lewisian  PDF,  gc  =  0.75.  Hence,  for  the 
Lewisian  model  field,  any  Gaussian  value  less  than 
0.75  is  mapped  to  gneiss,  and  everything  else  to 
amphibolite. 

2.  The  covariance  for  the  normal  equivalent  field, 
Cgg(x),  is  determined  such  that  when  the  PDF 
mapping  described  above  is  performed,  the  model 
field  conforms  to  the  model  covariance  Q,/,(x). 
Where  both  normal  equivalent  and  model  fields  have 
zero  mean  and  unit  variance,  the  relationship  is  given 
by  Christakos  [1992;  page  332,  equation  3],  In 
general  this  equation  must  be  solved  by  numerical 
integration. 

3.  Generate  a  synthetic  normal  equivalent 
realization  g(ij)  conforming  to  covariance  Cgg(k,l ) 
by  Fourier  methods  [e.g.,  Goff  and  Jordan,  1988; 
1989a]. 

4.  For  each  discrete  value  of  the  normal  equivalent 
g(i,j),  determine  the  value  of  the  model  field  h(ij) 
using  the  PDF  mapping  described  in  step  1 . 

The  Lewisian  Data/Synthetic  Comparison, 
Take  1 

The  Lewisian  2-D  stochastic  model  is  summarized 
as  follows:  (1)  The  PDF  is  a  binary  distribution  of 
78%  gneiss  (6.2  km/s)  and  22%  amphibolite  (6.75 
km/s);  (2)  the  covariance  is  modeled  by  a  von 
Kdrmdn  function  specified  by  parameters  kn  =  24.5 
km'  ,  ks  =  6.0  km"",  6S  =  85.0°  (where  0°  represents 
vertical  lineament  orientation),  and  D  =  2.54.  Figure 
lb  displays  a  full  two-dimensional  synthetic 
realization  of  the  stochastic  model  estimated  from  the 


Lewisian  map  of  Figure  la.  The  covariance 
parameters  for  the  normal  equivalent  field  are:  kn  = 
21.6  km'1,  ks=5.28  km'1,  05=  85.0°,  and  D  =  2.00. 

As  stated  in  the  Introduction,  comparison  between 
Figures  lb,  the  synthetic  realization  of  the  stochastic 
model,  and  Figure  la,  the  data  from  which  the 
stochastic  model  is  derived,  indicates  both  successes 
and  failures  of  the  stochastic  model.  The  success  of 
the  model  is  that  it  correctly  describes  the  physical 
properties  that  it  is  designed  to  describe;  i.e., 
probability  of  units,  scales  of  structures,  orientation, 
and  roughness.  The  failure  of  the  model  is  that  it  is 
not  complete;  i.e.,  it  cannot  reproduce  the  pattern 
which  is  visually  obvious  to  the  eye:  the  sinuous 
connectivity  of  the  units. 

III.  Sinuous  Connectivity  Modeling  and 
Synthetics 

Till  now  we  have  employed  a  stochastic  modeling 
technique  which  requires  joint  characterization  of  the 
PDF  and  covariance  of  the  field  of  interest.  A  great 
advantage  of  this  technique  is  that  estimation  of 
model  parameters  can  be  performed  through  direct 
inversion  of  sample  data.  An  entirely  equivalent, 
though  indirect,  stochastic  modeling  technique  would 
be  joint  specification  of  a  PDF  mapping  and  the 
covariance  function  of  a  normal  equivalent:  i.e.,  the 
"recipe"  for  generating  synthetic  realizations.  In  the 
case  of  the  Lewisian  model  discussed  above,  this  type 
of  stochastic  model  would  be  specified  by  the  normal 
equivalent  covariance  parameters  kn  =  21.6  km"1 ,  ks 
=  5.28  km"1,  6S  =  85.0°,  and  D  =  2.00,  and  the 
Gaussian-to-binary  PDF  mapping  function  with 
separation  scale  gc  =  0.75. 

The  rationale  for  adopting  the  latter  approach  to 
stochastic  modeling  is  that  it  provides  an  added 
dimension  to  characterization  not  available  in  the 
former:  the  PDF  mapping.  The  PDF  mapping 
presented  earlier  is,  in  fact,  only  the  simplest  possible 
case  of  a  Gaussian-to-binary  mapping.  The  only  real 
constraint  that  is  placed  on  this  mapping  is  that  the 
total  probability  of  unit  1  (gneiss)  is  0.78,  and  the 
total  probability  of  unit  2  (amphibolite)  is  0.22. 
Instead  of  finding  a  single  value  gc  which  separates 
the  Gaussian  PDF  into  regions  of  0.78  and  0.22 
probability,  we  can  formulate  a  more  complex 
mapping  into  alternating  bands  of  unit  1  and  unit  2, 
ensuring  that  the  sum  of  all  unit  1  probability  is  0.78, 
and  the  sum  of  all  unit  2  probability  is  0.22.  This 
operation  will  map  the  alternating  units  onto  contour 
sets  of  the  normal  equivalent  surface;  we  therefore 
apply  the  term  "contour  set  mapping"  to  this  type  of 
PDF  mapping.  As  anyone  who  has  seen  a  contour 
map  of  topography  will  quickly  intuit,  contour  set 
mapping  should  generate  a  binary  surface  possessing 
the  property  of  sinuous  connectivity. 

Contour  Set  Mapping  Synthetics 

More  formally,  we  specify  a  binary  contour  set 
mapping  in  the  following  way: 
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#2*  53  <  ^<xl>  <  ^4 

tfj,  52  <5(Xl)<  53 

#2-  5l  <  5(*1 )  <  52 

H{,  5o  =  “°°  <  5(*i)  <  5i 


(5) 


where  H i  and  H'l  represent  the  values  given  to  units 
1  and  2  respectively.  Each  interval  defined  by  the 
boundary  values  g(.  j  to  gj  is  defined  as  a  "contour 
set",  and  the  difference  g{  -  gi-i  is  defined  as  the 
"contour  set  thickness”.  Total  probabilities  for  each 
unit  are  computed  by  the  following: 


ft 

Pi=^  \pg(8'W> 


p2 = i  -  pi = y1,  j  Pcig'w 

ieven 

The  boundary  values  gi  must  be  chosen  so  that  pi 
and  p2  satisfy  our  total  probability  constraints  (for 
the  Lewisian  case,  pi  must  equal  0.78  and  p2  must 
equal  0.22).  Beyond  that,  however,  the  values  of  gi 
are  arbitrary,  unfortunately  resulting  potentially  in  an 
infinite  number  of  additional  parameters  to  the 
model.  Limitations  on  parameters  will  be  discussed 
below. 

While  contour  set  mapping  must  satisfy  total 
probability  constraints,  the  covariance  of  the  normal 
equivalent  field  must  be  chosen  so  that  the  model 
field  corresponds  to  the  model  covariance.  While  we 
cannot  at  this  time  make  formal  statements  of 
uniqueness,  it  is  nevertheless  clear  from  our 
experience  in  forward  modeling  that,  for  a  given  PDF 
mapping,  this  is  a  strong  constraint. 

Before  presenting  examples  of  synthetic 
realizations  generated  with  contour  set  mappings,  we 
must  remark  that  von  Kdrmdn  normal  equivalent 
covariances  do  not  strictly  imply  von  Klrmdn  model 
covariances,  as  was  the  case  for  simple  Gaussian-to- 
binary  mappings  [Goff  et  al.,  1994],  This  is 
particularly  true  where  contour  sets  for  either  unit  1 
or  unit  2  or  both  are  not  of  uniform  thickness.  Where 
unit  1  thickness  are  uniform  and  unit  2  thicknesses 
are  uniform,  the  model  covariance  resulting  from  a 
von  Kdrmdn  normal  equivalent  covariance  can  be 
reasonably  approximated  by  a  von  Kdrmdn  function. 
The  correspondence  between  the  model  covariance 
and  a  best-fitting  von  Kdrmin  function  in  these  cases 
is  not  always  exact,  but  misfit,  where  it  occurs,  is 
concentrated  at  the  larger  lags  where  resolution  of  the 
covariance  estimated  from  data  is  very  poor  (i.e„ 
either  will  fit  the  data  just  as  well).  However,  in 
some  cases,  two  of  which  will  be  highlighted  in  the 
following  section,  the  covariance  estimated  from  the 
data  field  is  not  well-fit  by  the  von  Kdrmdn  function; 


instead  of  a  simple  decay  with  increasing  lag,  the 
covariance  of  these  fields  exhibits  two  scales  of 
decay.  The  flexibility  allowed  in  the  contour  set 
mapping  is  instrumental  in  reproducing  these  more 
complicated  covariance  structures. 

Figure  3  displays  a  series  of  four  model  fields 
generated  using  progressively  thinner  contour  sets, 
labeled  as  model  fields  1-4.  The  contour  set 
mappings  are  also  shown,  displayed  graphically  as 
black  and  white  areas  below  the  Gaussian  PDF.  Each 
model  field  PDF  and  covariance  correspond  to  the 
stochastic  model  described  earlier  for  the  Lewisian 
map  (Figure  la).  Model  field  1  is  just  a  simple 
Gaussian-to-binary  mapping  presented  earlier  as 
Figure  lb.  All  the  normal  equivalent  fields  were 
generated  with  identical  random  phase  spectra,  which 
enhances  visual  recognition  of  structural  similarities. 

Our  primary  observation  derived  from  Figure  3  is 
that  contour  set  mappings  successfully  produce 
sinuous  connectivity  in  the  model  field.  In  particular, 
model  field  3  is  identical  to  the  synthetic  realization 
in  Figure  lc,  which  was  compared  to  the  data  field 
(Figure  la)  in  the  Introduction.  As  noted  earlier,  this 
comparison  is  visually  superior  to  the  comparison  of 
Figure  lb  (simple  Gaussian-to-binary  mapping)  to 
Figure  la. 

The  Sinuous  Connectivity  Parameter 

Now  that  we've  established  a  method  for 
characterizing  and  synthesizing  fields  that  possess  the 
property  of  sinuous  connectivity,  we  ask  ourselves  if 
it  is  possible  to  quantify  sinuous  connectivity  itself. 
In  other  words,  might  it  possible  to  state  that  one  field 
has  more  or  less  sinuous  connectivity  than  another 
field?  By  stepping  through  progressively  thinner 
contour  set  thicknesses.  Figure  8  was  designed  to 
address  this  question.  If  sinuous  connectivity  were  a 
quantifiable  property,  then  we  would  expect  a  general 
increase  in  this  property  as  we  decrease  in  contour  set 
thickness  in  the  PDF  mapping.  However,  this  does 
not  appear  to  be  the  case.  While  there  is  a  strong 
contrast  in  character  going  from  the  model  field  1  to 
model  field  2  in  Figure  3,  beyond  that  there  is  little 
noticeable  change  in  character  with  decreasing 
contour  set  thickness.  The  apparent  invariance  of 
sinuous  connectivity  characteristics  is  likely 
attibutable  to  the  self-similarity  property  of  the  fractal 
normal  equivalent  fields  [e.g.,  Mandelbrot,  1981]. 

The  important  implication  of  the  above 
observation  is  that  wide  classes  of  stochastic  models 
based  on  contour  set  mappings  and  normal  equivalent 
covariances  are  redundant,  or  non-unique.  On  the 
one  hand  this  is  problematic,  since  it  implies  that  this 
technique  for  stochastic  modeling  will  not  be  useful 
in  an  inversion  problem.  On  the  other  hand,  this 
observation  is  very  powerful,  because  it  vastly 
simplifies  the  parameterization  that  is  necessary:  the 
stochastic  model  is  complete  just  by  specifying  its 
PDF,  covariance,  and  whether  or  not  it  is  sinuously 
connective. 
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Model  Field  1 

Normal  Equivalent  Parameters: 

kn  =  21.6  km-1,  ks  =  5.28  km1,  D  =  2.00 


Model  Field  2 

Normal  Equivalent  Parameters: 

k„  =  4.90  km-1 ,  ks  =  1 .20  km'1 ,  D  =  2.04 


Model  Field  3 

Normal  Equivalent  Parameters: 

k„  =  1.89  km1,  ks  =  0.46  km1,  D  =  2.10 


Model  Field  4 

Normal  Equivalent  Parameters: 

kn  =  0.82  km-1,  ^  =  0.20  km1,  D  =  2.10 


Normalized  Gaussian  Variate 


Figure  3.  A  series  of  4  model  fields,  all  with  identical  PDF  and  covariance  (the  Lewisian  stochastic  model),  generated  with 
contour  set  mappings  of  decreasing  contour  set  thickness.  To  the  right  of  each  model  field,  normal  equivalent  covariance 
parameters  are.  All  normal  equivalent  fields  had  identical  H  (1.0)  and  9s  (85°).  Also  shown  are  to  the  right  are  graphical 
representations  of  the  contour  set  mapping;  any  value  sampled  from  the  normal  equivalent  field  that  fell  within  the  white 
regions  was  mapped  to  unit  1,  and  any  sample  that  fell  within  the  black  regions  was  mapped  to  unit  2.  In  Model  fields  2 
through  4,  the  contour  set  thickness  changes  by  a  factor  of  2  at  each  step.  Model  fields  1  and  3  are  identical  to  model  field 
presented  in  Figures  la  and  lb  respectively. 
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Figure  4.  (a)  Digitized  map  of  the  Franciscan  complex,  northern  California  (B.M.  Page,  1994,  unpublished).  A 
small  percentage  of  volcanic  units  within  the  map  were  ignored  as  insubstantial.  Map  is  66%  sandstone  and  34% 
melange.  While  sandstone  will  have  a  nearly  uniform  seismic  velocity  of  -5. 8-5.9  km/s,  melange  is  expected  to 
exhibit  a  range  of  velocities  from  -5.0-6. 4  km/s.  There  are  two  problems  in  stochastic  modeling  associated  with  the 
Franciscan:  the  interaction  between  melange  and  sandstone,  and  the  properties  of  the  melange.  For  the  present  we 
concentrate  on  the  former,  and  simply  identify  sandstone  and  melange  with  proxy  values  -1  and  +1.  (b)  Synthetic 
realization  based  on  PDF  and  covariance  modeling  of  the  Franciscan  complex,  including  2-scaled  contour  set 
mapping  and  the  property  of  sinuous  connectivity.  Normal  equivalent  parameters  are  listed  in  text,  (c)  Graphical 
representation  of  the  Franciscan  contour  set  mapping;  any  value  sampled  from  the  normal  equivalent  field  that  fell 
within  the  white  regions  was  mapped  to  unit  1  (sandstone),  and  any  sample  that  fell  within  the  black  regions  was 
mapped  to  unit  2  (melange),  (d)  Comparison  of  the  Franciscan  1-D  column-direction  (Chh(0j))  data  covariance 
(solid)  with  the  model  covariance  (dashed)  computed  for  the  2-scaled,  sinuously  connective  field. 
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Figure  5.  (a)  Digitized  map  of  the  Hafafit  complex,  eastern  Egypt  [Greiling  and  El-Ramly ,  1990;  Rashwan ,  1991]. 
The  gneissic  psammites  and  gneissic  granodiorites  have  only  incidental  contact,  and  so  were  considered  as  one  unit 
to  increase  coverage.  A  small  percentage  of  granites  were  ignored  as  insubstantial.  Where  it  was  clearly  obvious  to 
do  so,  units  were  interpolated  where  covered  by  wadi  alluvium.  Map  contains  50%  gneiss  and  50%  metagabbro. 
Petrophysical  data  are  not  known,  so  proxy  values  of  -1  and  +1  were  assigned  respectively  to  the  gneiss  and 
metagabbro.  It  is  certain,  however,  that  the  metagabbro  will  have  a  faster  seismic  velocity  than  the  gneiss,  (b) 
Synthetic  realization  based  on  PDF  and  covariance  modeling  of  the  Hafafit  complex,  including  2-scaIed  contour  set 
mapping  and  the  property  of  sinuous  connectivity.  Normal  equivalent  parameters  are  listed  in  text,  (c)  Graphical 
representation  of  the  Hafafit  contour  set  mapping;  any  value  sampled  from  the  normal  equivalent  field  that  fell 
within  the  white  regions  was  mapped  to  unit  1  (gneiss),  and  any  sample  that  fell  within  the  black  regions  was 
mapped  to  unit  2  (metagabbro).  (d)  Comparison  of  the  Hafafit  1-D  row-direction  (Qj/,(i,0))  data  covariance  (solid) 
with  the  model  covariance  (dashed)  computed  for  the  2-scaIed,  sinuously  connective  field. 
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IV.  Complex  Examples 

Two  more  examples  of  sinuous  connectivity 
modeling  are  presented  in  Figures  4  (Franciscan 
formation)  and  5  (Hafafit  gneiss  complex).  These 
examples  are  more  complex  than  the  Lewisian 
example  because  their  covariance  functions  were  not 
well-fit  by  a  von  Karman  model.  In  particular,  we 
required  two  superposed  scales  of  decay.  This 
behavior  can  be  matched,  however,  in  the  contour  set 
mapping  not  by  formulating  a  complex  normal 
equivalent  field,  but  rather  by  using  complex 
combinations  of  contour  thicknesses.  Normal 
equivalent  field  parameters  for  the  Franciscan  are:  kn 
-  0.144  km'  ,  ks  =  0.052  km"  ,  6S  =  -26.3°,  and  D  = 
2.00.  Normal  equivalent  field  parameters  for  the 
Hafafit  are:  kn  =  0.43  km'^,  ks  =  0.080  km  \  9S  = 
91.1°,  and  D  =  2.10. 

V.  Conclusions 

Previous  work  in  stochastic  modeling  of  lithologic 
heterogeneity  has  involved  joint  characterization  of 
the  PDF,  which  is  usually  modal  (i.e.,  binary,  ternary, 
etc.),  and  the  covariance  function,  which  is  often 
well-represented  by  the  von  Karman  function.  In  this 
paper  we  have  built  upon  that  earlier  work  by 
incorporating,  in  addition  to  PDF  and  covariance 
characterization,  a  method  for  modeling  fields  that 
possess  the  property  of  sinuous  connectivity.  This 
method  involves  defining  a  Gaussian-to-binary 
contour  set  mapping  and  a  normal  equivalent  field 
such  that,  when  the  normal  equivalent  field  is  mapped 
to  the  binary  model  field,  the  PDF  and  covariance  of 
the  model  field  honor  the  stochastic  model.  This 
modeling  scheme  also  constitutes  a  recipe  for 
generating  sinuously  connective  synthetic  realizations 
of  the  stochastic  model. 

An  important  observation  derived  from  synthetic 
realizations  is  that,  owing  to  the  self  similarity 
property  of  the  fractal  normal  equivalent  fields, 
sinuous  connectivity  modeling  is  highly  non-unique. 
Wide  classes  of  contour  set  mappings/normal 
equivalent  fields  will  generate  statistically  identical 
model  fields.  It  is  suggested,  therefore,  that  the 
property  sinuous  connectivity  can  be  characterized  by 
a  binary  parameter:  either  the  field  is  sinuously 
connective  or  it  is  not. 

Stochastic  models  and  synthetic  realizations  of 
lithologic  heterogeneity  should  play  a  critical  role  in 
modeling  seismic  wave  propagation  through  a 
complex  crust.  The  stochastic  model  presented  here 
is  the  most  realistic  and  complete  representation  of 
lithologic  heterogeneity  presented  to  date.  We 
believe  that  this  model  is,  in  fact,  sufficiently 
"complete"  for  application  to  the  seismic  problem; 
i.e.,  additional  complexity  or  parameterization  will 
probably  not  siginficantly  affect  the  observed  seismic 
wavefield.  Our  future  plans  center  around 
investigating  the  response  of  the  observed  wavefield 
to  parameters  of  the  stochastic  model. 
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